%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mapping reduced form coefficients into structural coefficients
% Identification of skill-biased and unskill-biased technology shocks using
% long-run sign restrictions
% this version: November 2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% inputs:
% BB                reduced form coefficients for ndraws
% Omega             VCV from reduced form residuals for ndraws
% nlags             number of lags
% ndraws            number of draws
% HORIZON_ident     identification horizon
% c                 choice of constant in the reduced form VAR
% num_vars          number of variables in the VAR
% nshocks           number of identified shocks
% hndx              index of hours in the VAR
% do_level          hours in levels of first differences

% output:
% Response          Responses to one standard deviation of all shocks for all horizons and all draws
% Atilde            matrix A from the structural representation
% Respres           Response to (all) one standard deviation residual shock(s)
%                   size is num_vars,HORIZON_ident,ndraws

function [Response,Respres]=ident_sign(BB,Omega,nlags,ndraws,HORIZON_ident,c,num_vars,nshocks,hndx,do_level)

Respdraws = [];
Respresdraws = [];
Adraws = [];
    
for n = 1:ndraws
    
    % Iteration
    dnum=num2str(n);
    if dnum(end)==num2str(0) && dnum(end-1)==num2str(0)
        disp('iteration'),disp(n)
    end
    
    % Step I:
    % transform the reduced coefficients into structural coefficents Phi
    Btilde = zeros(num_vars,num_vars,HORIZON_ident);
    for k=1:nlags
        Btilde(:,:,k) = BB(:,((k-1)*num_vars+c+1):(k*num_vars+c),n);
    end;
    Phi = zeros(num_vars,num_vars,HORIZON_ident+1);
    Phi(:,:,1) = eye(num_vars,num_vars); % this is time zero, i.e. the impact response
    
    for s = 2 : HORIZON_ident,
        Phitemp = zeros(num_vars,num_vars);
        for j= 1:(s-1),
            Phitemp = Phitemp + Phi(:,:,(s-j))*Btilde(:,:,j);
        end;
        Phi(:,:,s) = Phitemp;
    end;
    
    % Step II: Implementing the restrictions
    % i.e. get CCtilde and A, the identification horizon is flexible here
    % and may be adjusted
    CCtilde = eye(num_vars,num_vars);
    for j= 2:HORIZON_ident+1
        CCtilde = CCtilde + Phi(:,:,j);
    end;
    
    % we are now looking for some matrix A, such that AA'=Omega and such that
    % CCtilde*A is lower triangular
    Sigma=CCtilde*Omega(:,:,n)*CCtilde';
    
    % We will be decomposing the upper left of the long-run variance
    Slong_sign = Sigma(1:2,1:2);
    % starting value for the decomposition and rotation can be any matrix
    % that fulfills orthogonality assumption, usually Cholesky factor of the VCV is taken...
    T=chol(Slong_sign)';
    
    maxdraws = 100; % number of draws for the rotations
    % grid
    gridspace = 2*pi/(maxdraws-1);
    grid=zeros(maxdraws,1);
    for k=1:(maxdraws-1)
        grid(k+1,1) = grid(k,1)+gridspace;
    end
    
    % starting the rotations
    for dr = 1:maxdraws
        sep = 0;
        % rotation matrix
        theta = grid(dr);
        Q = [ cos(theta), -sin(theta);
            sin(theta), cos(theta) ];
        % rotated matrix of long-run effects
        L_sign = T*Q;
        
        % check sign restrictions
        % first shock is the SBT shock
        % check first column first
        sbtshock = 1;
        ubtshock = 2;
        z1 = [L_sign(1,1)>0 L_sign(2,1)>0];
        zz1=all(z1); % true if all elements in z1 are nonzero
        if zz1==0
            mz1 = [L_sign(1,2)>0 L_sign(2,2)>0];
            mzz1=all(mz1);
            if mzz1==1
                sbtshock = 2;
                ubtshock = 1;
            end
        end
        % second shock is a combination of neutral and UBT shock
        z2 = [L_sign(1,ubtshock)<=0 L_sign(2,ubtshock)>=0];
        zz2=all(z2); % true if all elements in z1 are nonzero
        if zz2 == 1
            if zz1 == 1
                sep = 1;
            elseif mzz1 == 1
                sep = 1;
            end
        end
        
        L_sign = T*Q;
        
        % calculating the remainder of L
        L = zeros(num_vars,num_vars);
        L(1:2,1:2) = L_sign;
        % calculating the remaining elements of the frist two columns
        for ii = 3:num_vars
            L(ii,1:2) = Sigma(1:2,ii)'/L_sign';
        end
        % decomposing the remainder with Cholesky (need to correct the long-run
        % VCV)
        sub_ndx = num_vars-2;
        CorrMat = zeros(sub_ndx,sub_ndx);
        for i=1:sub_ndx
            CorrMat(i,:) = L(3:end,1)'.*L(i+2,1) + L(3:end,2)'.*L(i+2,2);
        end
        Slong_sub = Sigma(3:num_vars,3:num_vars)-CorrMat;
        L_sub = chol(Slong_sub)';
        L(3:end,3:end) = L_sub;
        
        % calculating Atilde
        A_sign = CCtilde\L;
        
        % save the good draws
        if sep == 1
            
            Avec = reshape(A_sign,1,num_vars*num_vars);
            Adraws = [  Adraws
                Avec ];
            
            % Step III: IRFs for all shocks (separately)
            Resp = zeros(num_vars*num_vars,HORIZON_ident);
            for shock_counter = 1:num_vars,
                shock_vector = zeros(num_vars,1);
                if shock_counter == 1
                    shock_vector(sbtshock) = 1;
                elseif shock_counter == 2
                    shock_vector(ubtshock) = 1;
                else
                    shock_vector(shock_counter) = 1;
                end
                CCneutilde = zeros(num_vars,num_vars);
                for time_counter = 1 : HORIZON_ident,
                    CCneutilde = CCneutilde + Phi(:,:,time_counter);
                    Resphelp1 = CCneutilde*A_sign*shock_vector; % for differences
                    Resphelp2 = Phi(:,:,time_counter)*A_sign*shock_vector; % for levels
                    Resp((shock_counter-1)*num_vars+1:(shock_counter-1)*num_vars+num_vars,time_counter) = Resphelp1;
                    if do_level==1
                        Resp((shock_counter-1)*num_vars+hndx,time_counter) = Resphelp2(hndx,:)*100;
                    end
                end;
            end;
            
            Respdraws = [  Respdraws
                Resp ];
            
            % IRFs for all residual shock (bunched together)
            Resresp = zeros(num_vars,HORIZON_ident);
            shock_vector = zeros(num_vars,1);
            shock_vector(nshocks+1:num_vars) = 1;
            CCneutilde = zeros(num_vars,num_vars);
            for time_counter = 1 : HORIZON_ident,
                CCneutilde = CCneutilde + Phi(:,:,time_counter);
                Resphelp1 = CCneutilde*A_sign*shock_vector; % for differences
                Resphelp2 = Phi(:,:,time_counter)*A_sign*shock_vector; % for levels
                Resresp(:,time_counter) = Resphelp1;
                if do_level==1
                    Resresp(hndx,time_counter) = Resphelp2(hndx,:)*100;
                end
            end;
            
            Respresdraws = [  Respresdraws
                Resresp ];

        end
    end 
    
end; % draws loop

gdraws = size(Adraws,1);
Atilde = zeros(num_vars,num_vars,gdraws);
Rsize = num_vars*num_vars;
Response = zeros(Rsize,HORIZON_ident,gdraws);
Respres = zeros(num_vars,HORIZON_ident,gdraws);
for g=1:gdraws
    % reshaping back A
    Atilde(:,:,g) = reshape(Adraws(g,:),num_vars,num_vars);
    % reshaping back the responses
    Response(:,:,g) = Respdraws((g-1)*Rsize+1:(g-1)*Rsize+Rsize,:);
    Respres(:,:,g) = Respresdraws((g-1)*num_vars+1:(g-1)*num_vars+num_vars,:);
end